Fusariose: QTLs Description - Part2
Fusariose: QTLs Description - Part2
This file follows the first part of the QTL analysis (Description_of_QTLs-Part1). Load the results + functions of the first part:
load("/Users/yan/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/QTL_R_environment.R")Charge some libraries that will be useful
library(RColorBrewer)
library(xtable)
library(plotly)
library(FactoMineR)7/ IC on the map (Fig)
We can also visualize the IC on the genetic map. Let’s visualize important variable ONLY = AUDPC and last note for BLUPs only.
par(mfrow=c(1,1))
my_var_to_show=PUBLI
show_IC_QTL_on_map(data, my_var_to_show , c("1B", "5A"), my_QTL_thres , 1.5)I save it as a png image for publication purpose.
png("~/Dropbox/Publi_Fusariose/FIGURE/Fig_IC_QTL_on_map.png")
show_IC_QTL_on_map(data, my_var_to_show , c("1B", "5A"), my_QTL_thres , 1.5)
dev.off()## quartz_off_screen
## 2
8/ Evolution of LOD scores with time?
We are going to study qtl 1B and 5A only. Let’s have a look to the maximum LOD score of these QTL at each date of observation. For each experiment, for each variable?
Functions
Fonction qui regroupe chaque LOD max dans un tableau. On ne garde que les BLUPs de toutes les variables à toutes les dates.
find_LOD_max=function( chromo, my_min, my_max ){
bilan=data.frame(matrix(0, 200 ,8))
num=0
vec_names=c("NOT", "PEPI", "NEPIL", "PEPIL")
for(expe in list(CAP13, GRI11, GRI13, GRI15, LEC14, BEQ11, BAR14)){
num_carac=0
for(carac in list(NOT, PEPI, NEPIL, PEPIL)){
num_carac=num_carac+1
my_carac=intersect(BLUP, intersect(expe, carac ))
if( expe==BEQ11){ my_carac=intersect(expe, carac ) }
if( length(my_carac)==0) { next }
num=num+1
my_vec=c()
for(i in c(1:length(my_carac))){
b=(data$LOD[ which(data$LG==chromo & data$variable==my_carac[i] & data$Distance<my_max & data$Distance>my_min ) ])
a=max( b , na.rm=T)
my_vec=c(my_vec, a)
}
my_replace=c( unique(gsub("_.*","",my_carac)) , vec_names[num_carac], my_vec )
bilan[num, c(1:length(my_replace)) ] = my_replace
}
}
bilan=bilan[bilan[,1]!=0 , ]
bilan[bilan==0]=NA
return(bilan)
}Fonction pour visualiser ce tableau.
show_max_LOD=function(bilan){
# Préparation des axes des X.
x1=c("300", "400", "550", "AUDPC")
x2=c("300", "350", "450", "500", "550", "AUDPC")
x3=c("350", "450", "550", "AUDPC")
x4=c("350", "550", "AUDPC")
x5=c("300", "400", "550", "AUDPC")
x6=c("300", "350", "450", "500", "AUDPC")
x7=c("300", "400", "550", "600", "AUDPC")
my_x_axis=list(x1,x2,x3,x4,x5,x6,x7)
par(mfrow=c(2,4) , mar=c(3,3,1,1))
liste_var=c("NOT", "PEPI","NEPIL", "PEPIL")
num=0
for(expe in list("CAP13", "GRI11", "GRI13", "GRI15", "LEC14", "BEQ11", "BAR14")){
num=num+1
AA=bilan[ which(bilan[,1]==expe),]
plot(seq(1,6) , AA[ 1, c(3:ncol(AA))] , type="l", ylim=c(0,10) , col=my_colors[match( AA[ 1, 2] , liste_var )], lwd=2, xlab=expe, ylab="", xaxt="n")
axis(1, at=seq( 1 , length(my_x_axis[[num]])), my_x_axis[[num]], cex=0.7)
abline(h=my_QTL_thres, col="grey")
abline(v=c(1:6), col="grey", lwd=0.3)
for(i in c(2:nrow(AA))){
points(seq(1,6) , AA[ i, c(3:ncol(AA))] , type="l" , col=my_colors[match( AA[ i, 2] , liste_var )], lwd=2)
}
text(x=1.8, y=9, expe, col="orange" )
}
plot(1,1,col="transparent", xaxt="n", yaxt="n", bty="n", xlab="", ylab="")
legend( "center", legend = c("NOT", "PEPI","NEPIL", "PEPIL") , bty="n", pch=20 , col=my_colors[1:4], pt.cex=4, cex=1.6 )
}QTL 1B
Calcul des LOD max à chaque date:
max_lod_1B=find_LOD_max( "1B", 0, 9 )Show it:
show_max_LOD(max_lod_1B)Donc ces graphiques représentent les valeurs des LODs max de chaque QTL, i.e. pour chaque variable de chaque expe pour le chromosome 1B.
QTL 5A
Calcul des LOD max à chaque date:
max_lod_5A=find_LOD_max( "5A", 100, 300 )Show it:
show_max_LOD(max_lod_5A)9/ Link mith physical reference
OK so we have a QTL with a confidence interval. It could be great to see how many genes are present in this IC, and the physical size it represents.
1B
Summary of the 1B IC of the QTL:
# initialise the table
bil_phy=data.frame(matrix(0,0,9))
colnames(bil_phy)=c("chromo","start IC (cM)", "end IC (cM)", "len IC (cM)", "# SNP", "start IC (Mb)", "end IC (Mb)", "size IC (Mb)", "# genes")
# 1B
bil_phy[1,]=find_IC_infos(1.2,14.1,"1B",c())| chromo | start IC (cM) | end IC (cM) | len IC (cM) | # SNP | start IC (Mb) | end IC (Mb) | size IC (Mb) | # genes |
|---|---|---|---|---|---|---|---|---|
| 1B | 1.2 | 14.1 | 12.5 | 58 | 1527986 | 6092666 | 4564680 | 0 |
This IC has 103 genes. Let’s check the colinearity between genetic and physical positions in this area.
#1B
par(mfrow=c(1,2), mar=c(5,5,3,3))
show_genes_of_IC(0, 300, "1B", 0.005, 1, c())
show_genes_of_IC(0, 14, "1B", 0.5, 4, c())5A
Summary of the 5A IC of the QTL:
#5A
bil_phy[2,]=find_IC_infos(235.6 ,305,"5A","Traes_5AL_99116D900@696")| chromo | start IC (cM) | end IC (cM) | len IC (cM) | # SNP | start IC (Mb) | end IC (Mb) | size IC (Mb) | # genes |
|---|---|---|---|---|---|---|---|---|
| 5A | 235.6 | 305 | 68.7 | 86 | 118853998 | 143979436 | 25125438 | 0 |
#5A
par(mfrow=c(1,2), mar=c(5,5,3,3))
show_genes_of_IC(0, 400, "5A", 0.005, 1, c())
show_genes_of_IC(235.6 ,305, "5A", 0.1, 4,"Traes_5AL_99116D900@696")4B (height)
Summary of the 4B IC of the QTL:
# 4B
bil_phy[3,]=find_IC_infos(47,53,"4B","")## Warning in min(IC_map$Distance): aucun argument trouvé pour min ; Inf est
## renvoyé
## Warning in max(IC_map$Distance): aucun argument pour max ; -Inf est renvoyé
## Warning in min(IC_map$Posi_physique): aucun argument trouvé pour min ; Inf
## est renvoyé
## Warning in max(IC_map$Posi_physique): aucun argument pour max ; -Inf est
## renvoyé
| chromo | start IC (cM) | end IC (cM) | len IC (cM) | # SNP | start IC (Mb) | end IC (Mb) | size IC (Mb) | # genes |
|---|---|---|---|---|---|---|---|---|
| 4B | 47 | 53 | -Inf | 0 | Inf | -Inf | -Inf | 0 |
#4B
par(mfrow=c(1,2), mar=c(5,5,3,3))
#show_genes_of_IC(0, 200, "4B", 0.005, 1)
#show_genes_of_IC(47, 53, "4B", 0.1, 4)7B (precocity)
Summary of the 7B IC of the QTL:
# 7B
bil_phy[4,]=find_IC_infos(5,35,"7B","")## Warning in min(IC_map$Distance): aucun argument trouvé pour min ; Inf est
## renvoyé
## Warning in max(IC_map$Distance): aucun argument pour max ; -Inf est renvoyé
## Warning in min(IC_map$Posi_physique): aucun argument trouvé pour min ; Inf
## est renvoyé
## Warning in max(IC_map$Posi_physique): aucun argument pour max ; -Inf est
## renvoyé
| chromo | start IC (cM) | end IC (cM) | len IC (cM) | # SNP | start IC (Mb) | end IC (Mb) | size IC (Mb) | # genes |
|---|---|---|---|---|---|---|---|---|
| 7B | 5 | 35 | -Inf | 0 | Inf | -Inf | -Inf | 0 |
#7B
#par(mfrow=c(1,2), mar=c(5,5,3,3))
#show_genes_of_IC(0, 200, "7B", 0.005, 1)
#show_genes_of_IC(5, 35, "7B", 0.1, 4)summary
| chromo | start IC (cM) | end IC (cM) | len IC (cM) | # SNP | start IC (Mb) | end IC (Mb) | size IC (Mb) | # genes |
|---|---|---|---|---|---|---|---|---|
| 1B | 1.2 | 14.1 | 12.5 | 58 | 1527986 | 6092666 | 4564680 | 0 |
| 5A | 235.6 | 305 | 68.7 | 86 | 118853998 | 143979436 | 25125438 | 0 |
| 4B | 47 | 53 | -Inf | 0 | Inf | -Inf | -Inf | 0 |
| 7B | 5 | 35 | -Inf | 0 | Inf | -Inf | -Inf | 0 |
10/ Comparaison avec les QTLs de la biblio
Using SSR
Let’s have a look to all the 9 SSR markers I have in my map:
Is it normal I have only 9?
This is not gonnna help a lot…
- 1B: pas de SSR
- 5A: 2 SSR, mais très très loin de l’IC.
–> Il faut aller voir si je peux trouver des arqueurs de la littérature qui ont été passé sur Dic2 x Silur, et les replacer sur la carte comme pour la publi mosaiquE…
SSR=unique(data[ grep("gwm" , data$marqueur) , c(1:3) ])
nrow(SSR)## [1] 9
| LG | marqueur | Distance |
|---|---|---|
| 1A | Xgwm164 | 68.40 |
| 3B | Xgwm376 | 92.30 |
| 3B | Xgwm77 | 94.20 |
| 4B | Xgwm495 | 94.00 |
| 5A | Xgwm205 | 41.00 |
| 5A | Xgwm415 | 72.50 |
| 5B | Xgwm544 | 55.30 |
| 7A | Xgwm276 | 181.50 |
| 7B | Xgwm400 | 38.50 |
Position of Snn1 on map
Le gène Snn1 est positionné sur le 1B et apporte une résistance forte à la septoriose, un autre champignon. On peut trouver la séquence de ce gène ici:
https://www.ncbi.nlm.nih.gov/nuccore/858936371/ Publié par: http://advances.sciencemag.org/content/2/10/e1600822.full
Blast sur la référence blé tendre release 28
makeblastdb -in sequence_snn1.fasta -dbtype nucl -out database
Blast à 100% sur ce contig: Traes_1BS_C4B727779
Ce contig n’est pas dans notre carte génétique, donc on a pas de SNP dessus. Par contre, en terme de position physique, il se situe à 880274 Mb sur le 1BS. Ca correspondrait donc à être entre 1.19 et 1.59 cM sur notre carte.
1B Traes_1BS_7BC6EFB55@1326 1.19999999999999 1B 502322 1B Traes_1BS_81F7DAE3E@432 1.59999999999999 1B 2538668
Nous notre QTL a tendance a etre un poil plus loin, mais vraiment vraiment proche…
11/ Relationship QTL presence / field characteristic
The intensity of the Fusarium infection was not the same from one experiment to another. Let’s check the relationship between disease intensity and strength of the QTL detected?
12/ Epistatic effect?
OK so we have 2 main QTLs. Do we have a epistatic interaction betwwen them??
Example
Let’s start on an example: BAR14_PEPIL300_blup. 2 QTLs: markers Cluster_16778|Contig1|original@1840 et Cluster_9940|Contig1|complementarySeq@104
# calculate inter
bil_inter=data.frame(matrix(0,0,9))
colnames(bil_inter)=c("trait","marker-1B", "marker 5A", "mean A-A", "mean A-B", "mean B-A", "mean B-B", "R2 tot", "pval inter")
bil_inter[1,]=analyse_inter("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup")| trait | marker-1B | marker 5A | mean A-A | mean A-B | mean B-A | mean B-B | R2 tot | pval inter |
|---|---|---|---|---|---|---|---|---|
| BAR14_PEPIL300_blup | Cluster_16778|Contig1|original@1840 | Cluster_9940|Contig1|complementarySeq@104 | -0.2 | 3.32 | -0.91 | -0.65 | 0.281 | 0.0001369779 |
And show this interaction!
boxplot_two_QTL_interactive("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup")#boxplot_two_QTL("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup", ylim=c(-4,7))All trait
Il faut commencer par trouver tous les traits qui ont un QTL signif sur le 1B et un signif sur le 5A.
var_with_2QTLs=c()
for(i in list(summary_QTL_NOT, summary_QTL_PEPI, summary_QTL_NEPIL, summary_QTL_PEPIL)){
a=i[which(i$chromo%in%c("1B","5A")) , ]
b=table(a$carac)
var=names(b[b==2])
var_with_2QTLs=c(var_with_2QTLs, var)
}
var_with_2QTLs=var_with_2QTLs[grep("blup",var_with_2QTLs)]| var_with_2QTLs |
|---|
| BAR14_PEPI350_blup |
| BAR14_PEPI400_blup |
| BAR14_PEPI600_blup |
| BAR14_PEPIAUDPC_blup |
| GRI13_PEPI550_blup |
| BAR14_NEPIL300_blup |
| BAR14_PEPIL300_blup |
Calculons les p-values de l’effet d’interaction pour tous ces cas:
for(i in var_with_2QTLs[-which(var_with_2QTLs=="BAR14_PEPIL300_blup")]){
vv=analyse_inter("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , i)
bil_inter=rbind(bil_inter , vv)
}On a donc un effet d’épistasie qui est significatif à tous les coups (quasiment).
Représentons ces épistasies avec un boxplot pour chacune des variables:
# widget to choose data
#selectInput("select", label = h3("Select a variable"), choices = var_with_2QTLs, selected = 1)
# interactibe boxplot
#renderPlotly({
# boxplot_two_QTL("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , select)
#})Dans la majorité des cas, on observe que les individus avec 1 seul QTL ou avec 2 QTLs sont similaires. Par contre, les individus avec aucun des 2 sont très différents et très sensibles. Donc en gros, si j’ai au moins un des 2 gènes tout va bien, sinon je suis sensible. Donc on est pas dans un modèle additif à priori, les 2 gènes sont liés.